1/u;- flux-noise and dynamical critical properties of two-dimensional XY-models 
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We have numerically studied the dynamic correlation functions in thermodynamic equilibrium of 
two-dimensional 0(2)-symmetry models with either bond (RSJ) or site (TDGL) dissipation as a 
function of temperature T. We find that above the critical temperature the frequency dependent 
flux noise ~ |l + (w /fi) 2 |~ Q(T)/2 , with 0.85 < a(TDGL)(T) < 0.95 and 1.17 < a(RSJ)(T) < 

1.27, while the dynamic critical exponents z(TDGL) ~ 2.0 and z(RSJ) ~ 0.9. Contrary to 
expectation the TDGL results are in closer agreement with the experiments in Josephson-junction 
arrays by Shaw et at, than those from the RSJ model. We find that these results are related to 
anomalous vortex diffusion through vortex clusters. 
PACS numbers: 74.20.-z, 64.60. Fr, 74.50. +r, 74.40. +k 



There is a general theoretical consensus that the equi- 
librium critical properties of two-dimensional Abelian 
symmetry statistical mechanical models are properly 
described by the Berezinskii-Kosterlitz-Thouless (BKT) 
theory [1,2]. The theory is based on the idea that at 
low temperatures there are thermally excited vortex- 
antivortex pairs that unbind just above the BKT critical 
temperature. Experimental tests have for the most part 
been indirect, in out of equilibrium conditions. Successful 
comparisons of the BKT predictions to experiment were 
initially carried out in superfluid helium [3] and in super- 
conducting films [4]. Later measurements on Josephson- 
junction arrays (JJA) yielded current-voltage character- 
istics [6], and dynamic impedance [7] results that also 
agree with the BKT theory. To explain the experiments, 
the static BKT scenario was heuristically extended to 
non-equilibrium situations [5] . One of the basic tenets of 
the dynamical BKT extensions is that vortices move dif- 
fusively with binding and unbinding of vortex pairs under 
the action of the external drive. Further phenomenolog- 
ical modifications to the dynamic extension of the BKT 
theory, that reduced the number of fitting parameters, 
have led to better correspondence to some of the experi- 
mental results on superconductors [8]. 

A very recent flux noise SQUID experiment by Shaw et 
al. [9] on proximity effect JJA provided a direct experi- 
mental test of the equilibrium BKT theory, at Tbkt < T 
[10]. They found the flux noise to be white at the lowest 
frequencies, and proportional to \/uj a , with a R* 1, at 
intermediate ones rather than the 1/lo 2 expected by the 
phenomenological theories [5,8]. They also carried out a 
dynamical scaling analysis of the noise function for dif- 
ferent temperatures, that yielded a dynamical exponent 
of z ~ 2. An earlier, less extensive, JJA experiment by 
Lerch et al. [11] found al/w behavior for the T < Tbkt 
regime, without the white noise frequency region. In this 



paper we concentrate on the Tbkt < T regime. We show 
that by numerically calculating quantities closely related 
to those measured in [9] we do get good agreement with 
their a and z results for a site dissipation model defined 
below [12]. 

There are two models that have been considered to 
study the arrays. One is the resistivcly-shuntcd Joseph- 
son (RSJ) array model and the other is related to the 
time-dependent Ginsburg-Landau (TDGL) model. The 
RSJ model is constructed from the elementary RSJ equa- 
tions for single Josephson junctions that form the array 
units, plus Kirchhoff's current conservation conditions 
at each lattice site. The JJA-RSJ model has been suc- 
cessful in explaining, for example, the experimental gi- 
ant Shapiro steps that arise when the JJA is driven by a 
dc+ac current [13]. The TDGL model is an alternative 
dynamical model that has been used to describe the gen- 
eral critical dynamical properties of the XY-model [14] 
and also some aspects of JJA behavior, including its flux 
noise [15,16]. An important theoretical question is then 
to ascertain which model is best to describe the exper- 
imental noise results. To find out, we have studied the 
flux noise and dynamical scaling properties for both mod- 
els. We do indeed find that the flux noise is anomalous 
in both cases, but the TDGL results appear to be closer 
to the experimental ones. We find that the vortices do 
not move independently but in a sea of vortex clusters 
that modify their diffusive properties. 

The equation of motion for both models can be written 
in the general Langevin equation form, 

f)M\0] 

d t d(r, t) = T(r, r')^^ + W(r, r', t) }. (1) 

Here 6(r,t) is the phase of the order parameter of the 
superconducting island at site r; H[6] = Ej cos(#(r', t) — 
0(r,t)) is the Hamiltonian bond energy, with Ej the 
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Josephson coupling, and W(r, t; r', t') the noise function. 
The RSJ model has T RS j{r,r') = rG(r,r'). In our 
units r = 1. G(r, r'), the two-dimensional inverse lat- 
tice Laplacian, arises from the fact that the currents in 
the array are conserved at each lattice site. The dis- 
sipation here is present in the junctions between the 
superconducting islands. The random noise function 
Wrsj{ v , t; r', t') = rj(r, t; r', t') is defined at each bond in 
the lattice with Gaussian properties < r)(r, t;r',t') >= 
and < r?(r, t; r + e i; t')r](r', t; r' + e 3 ,t') >= 2T5ij5(r - 
r')S(t — t'), where ej is a unit vector along the jth direc- 
tion, and T is the dimensionless temperature. 

The TDGL model is defined by taking T TDGL (r, r') = 
T8{r - r') (T = 1), and W TDGL (r, t; r', f ) = r,(r, t)5(r - 
r')S(t — t'); here the noise is also white and defined at 
each island, with < rj(r,t) >= and < 77(1* , t)r){r' , t') >= 
2T5(r — r')S(t — t'). This model has been used before to 
study the critical dynamics of the XY-model. This is in 
part because it is computationally less demanding than 
the RSJ model, that needs time consuming matrix inver- 
sions to evaluate G(r, r'). We find that the macroscopic 
thermodynamic properties, like the hclicity modulus and 
the phase correlation functions described below, are es- 
sentially the same for both models. However, as we show 
here, other dynamical equilibrium properties are differ- 
ent for the two models, and much more so the out-of- 
cquilibrium quantities [17]. 

An important aspect of the Shaw et al. study is that 
the SQUID used for measurements was smaller than the 
array size. The SQUID measures the average net num- 
ber of thermally nucleated and annihilated vortices that 
randomly enter or exit its effective area. We do exactly 
the same in our calculations, by considering regions of 
the lattice smaller in size than that of the full array sim- 
ulated (see Fig. 1(b)). The space-time local vorticity at 
plaquette R is defined by 27m(R, t) = J2v(R) [^( r ' — 
9(r', t)] mod(27r). Here T'(R) denotes an anti-clockwise 
sum around plaquette R. We imposed periodic boundary 
conditions (PBC) in both directions in our equilibrium 
calculations. The average vorticity for the ith SQUID of 
area I x I is defined as N\{£) — J2ne(ixi)i n (R>i)- To 
get the net flux noise we average over all the i-SQUIDS 
that can be fitted in the array of size L x L (See Fig. 
1(d).) The flux noise produced by the vortices in the 
time domain is then generally defined by 

= ]^Ei wmm - w(o)) 2 ]. (2) 

i 

Here Ni is the total number of i-SQUIDS of size £ x £, 
and the averages <,> are carried out over the prob- 
ability distribution for the noise or over time. Based 
on the ergodic theorem, both averaging methods should 
and do give the same results, but the latter is easier to 
implement numerically. The experimentally measured 
flux noise is simply given by the Fourier transform of 



gj{t) i.e. S*(w) = J dt e iult gj{t). To analyze the scal- 
ing results for 5$ we also need to calculate the time- 
dependent equilibrium phase correlation function, de- 
fined as g e {r,r',t,t') =< e W(^t)-e(r',t')] > In our ex _ 

plicit calculations we evaluate the zeroth-momentum cor- 
relation function, that is known to have only one dom- 
inant Lyapunov exponent [18]. Here we concentrate on 
the Tbkt < T region, as was done in the Shaw et al. 
experiment. In the long distance regime 

g e (r)=Ar- r >e- r ^ T \ (3) 

where the correlation length £ is expected to have the 
BKT form, 

£(T) = ( e V/J " TGKT , (4) 

with Tbkt, Co an d b constants determined in the cal- 
culations. In Fig.l we show the results for the equilib- 
rium time-averaged ge{r) for different temperatures, for 
both the TDGL (a) and RSJ (c) models in a 64 x 64 ar- 
ray. In the insets we plot the corresponding £(T) data 
obtained from the correlation function calculations, to- 
gether with their BKT-type fits. The fitting parameters 
are Co = 0.230(1), b = 1.45(1) and T BKT = 0.935(1) (in- 
set in Fig.l(a) TDGL) and Co = 0.274(30), b = 1.41(10) 
and Tbkt = 0.917(10) (inset in Fig. 1(c) RSJ). These 
non-universal results compare well with the numbers ob- 
tained from a more extensive Monte Carlo calculation 
that gave Co = 0.205, b = 1.6113 and Tbkt = 0.9035 
[18]. The reasonable agreement between these and the 
MC results gives strong support to the reliability of our 
dynamic simulation algorithms [19]. 

The results for gj (t) are shown in Fig. 1 (b) (TDGL) 
and 1 (d) (RSJ), for an i = 32 SQUID. As in the ex- 
periment we have restricted the analysis to the temper- 
ature regime with short correlation lengths in the range 
C < I < L. We tried different reasonable fitting func- 
tions for the data for gj (t) in the short and intermediate 
time regimes. From stretched exponentials, that gave 
reasonable results, to a novel and somewhat compara- 
ble fit, with fewer parameters, that involves the modified 
Bessel function K v (t), with v close to zero. The moti- 
vation for using the latter functions will become clearer 
below. First we show the noise functions S<s>(w) in Fig. 2 
(TDGL) (a) and (RSJ)(c). These figures were obtained 
from a direct fast Fourier transform of the gj{t) data. 
We see there that the intermediate frequency regime of 
5$ can be approximately fitted by a straight line lead- 
ing to 5$ ~ l/\oj\ at " T \ The a exponent has a weak 
non-monotonic temperature dependence but within a 
rather narrow range, 0.85 < a(TDGL){T) < 0.95, and 
1.17 < a(RSJ)(T) < 1.27, as shown in the inset to Figs. 
2 (a) and (c). At lower frequencies there is a bending of 
S<p tending towards a flat behavior for the lower temper- 
atures considered. We note that the TDGL a values are 
closer to the experimental values which were about one. 
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Although the RSJ results have slightly larger a's they 
can still be classified as anomalous l/\u>\ noise. Based on 
these results we propose an analytic ansatz for 5$(w) of 
the form S$(lj) ~ |£1 2 + cj 2 |~ q ( t )/ 2 . Here Q, is a charac- 
teristic frequency which in our calculations appears to be 
small. Further justification for this ansatz will be given 
below. We can then obtain the corresponding expression 
for gj (t) from the inverse Fourier transform of the ansatz 
for S<s>(u) giving, 

gY(t) = Ct» r ^l ) /2) K v (nt), (5) 

where v = (a(T) — 1)/2, T(x) is the Gamma function and 
C is a constant. We have used this result to fit the numer- 
ical data for gj and the fits are shown as continuous lines 
in Figs. 1(b) and 1(d). Note that based on the calculated 
values for a, the Bessel function index v takes the values, 
-0.075 < v(TDGL) < -0.025 and 0.085 < is(RSJ) < 
0.135 (recall that K_ v {x) = K v {x)). It is clear that the 
TDGL model is rather close to the Kg result, while the 
RSJ is slightly different. To further understand our nu- 
merical results, given that v is rather close to zero, it is 
instructive to first consider the case when a = 1, which 
leads to gj oc K (SU) oc e - nt J °° <It[t(t + t)]~ 1/2 , 
with Fourier transform oc (f2 2 + |kj| 2 )~i. In the limit 
when iv 3> f2, S<s,(uj) ~ l/|w|, while in the opposite limit 
it is a constant and the noise is white, just as in the 
experiment. The large u> regime corresponds to the log- 
arithmic short time region of K (Qt), whereas the lower 
uj regime corresponds to the exponential one. We can go 
further by looking at the definition of gj given in Eq.(2), 
which implies a coarse grained N{t) ~ i _1 ' 2 e _n * (for 
a = 1). Here we are making use of the ergodic theorem. 
We can extend the analysis given above to v ^ 0, which 
yields N(t, v{T)) ~ £-(V2-") e -fM. The evaluation of the 
time Fourier transform of < N(t, v{T))N(t + t, v{T)) > 
gets us back to our ansatz for S$(u>). This is an interest- 
ing result that can be interpreted in terms of an incoher- 
ent superposition of a series of independent random wave 
front events produced by the entering and leaving of vor- 
tices from the effective SQUID area. The coarse grained 
events can be assumed to be Poisson distributed in time 
so as to produce a random pulse train with the explicit 
time dependence for g v given above (see Ref. [20]). 

Another important experimental result was the test 
of the dynamic scaling predictions in frequency space. 
We have carried out an equivalent scaling analysis in 
time. The general dynamic scaling hypothesis assumes 
that gj = ^F^t/r^L /£,,£/£). We find that, as in the 
experiment, for £ < I < L, we can collapse the data into a 
single curve by considering only the time dependent part 
of the scaling function, i.e. g( cx £ /3 G(t/r^), with (3 an 
exponent and the relaxation time oc £ z , thereby defin- 
ing the dynamical critical exponent, z. We simulated a 
64 x 64 array for different temperatures and we were able 



to collapse all the data into one single curve for differ- 
ent T's and fs using this criterion. The data collapsed 
corresponds to roughly a decade in t£. We used data for 
7 different T's with £ = 16 and i = 32, for £(T)e[2.2, 8]. 
The accessible length scales and time scales are limited 
by the available computer power. For a given £, we used 
an L > 8£ lattice and simulated for 100-lOOOr^. The 
scaling results are shown in Fig. 2(b) (TDGL) and 2(d) 
(RSJ). From the scaling analysis we determined z and 
obtained the values of z(TDGL) - 2 and z(RSJ) - 0.9, 
as shown in the insets to Figs. 2(b) and 2(d). The es- 
timated z(TDGL) ~ 2 exponent is consistent with the 
value obtained in the experiment [9]. 

We have also made animations of the vortex motions 
above Tbkt- There we see that the dynamics is domi- 
nated by vortex clusters formed above Tbkt- One can 
thus link the anomalous values of the a exponents to the 
anomalous way in which single vortices diffuse through 
the clusters. 

In conclusion, we have studied the flux noise in the 
TDGL and RSJ models with the goal of providing an 
understanding to the recent interesting experiments by 
Shaw et al. [9]. We find that by following the prescrip- 
tions defined in the experiment, both models lead to 
anomalous vortex dynamics, with the TDGL results be- 
ing in closer agreement to the experimental ones. It is 
possible that the actual samples that are made of large 
area Nb islands deposited on a copper substrate have a 
very small resistance to ground and thus the TDGL com- 
ponent may dominate the bond or RSJ term [22] . We also 
have studied the effect of a small magnetic field and we 
find that the results change significantly. The a expo- 
nent acquires a stronger temperature dependence. Fur- 
thermore, we have preliminary results for the T < Tbkt 
regime and the a-values obtained are also indicative of 
anomalous vortex dynamics. We will report on these and 
other results in a separate publication. 
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FIG. 1. Results from simulations of a 64 x 64 array 
with PBC. Phase correlation functions ge{x) obtained from 
TDGL(a) and RSJ(c) simulations. The different curves cor- 
respond to: T = 1.08 (o), T = 1.12 (squares), T = 1.24 
(diamonds) and T = 1.32 (A). The continuous lines are fits 
to Eq. (3). Insets to (a) and (c) give the £(T) results for 
TDGL and RSJ dynamics, with the continuous lines fits to 
the BKT form given in Eq.(4), respectively. In (b) and (d) 
we plot gY (t) for £ = 32. The temperatures and notation is 
the same as in (a) and (c). The solid lines are fits to Eq.(5). 

FIG. 2. Flux-noise S<s>(lu) vs uj plotted on log- log scale for 
TDGL(a) and RSJ(c) dynamics for the same parameters as 
in Fig. 1. The different temperature curves have been shifted 
upwards by a constant distance for clarity. From top to bot- 
tom (a) T = 1.08, 1.16, 1.24, 1.32, 1.40 and 1.48 and (c) 
T = 1.08, 1.10, 1.12, 1.20,1.28 and 1.36. In the insets to pan- 
els (a) and (c) we plot the exponent a(T) obtained from a least 
squares fit to the quasi-linear frequency regime. The scaling 
plots for gY(t) are shown in panels TDGL(b) (T = 1.12, 1.16, 
1.20, 1.24, and 1.28), and RSJ(d) (T = 1.12, 1.16, 1.20, 1.24, 
1.28, 1.32, and 1.36), with £=16 (+) and £ = 32 (o). In the 
insets to panels (b) and (d) we show the scaling of t£ vs £ on 
a log-log plot, for £ — 32 (o) and I = 16 (square). The con- 
tinuous (dot-dashed) line is a linear fit to the £ — 32 (£ = 16) 
results. 
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